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1 Introduction 

Eddy current calculation is an important subject in classical electrody- 
namics. Eddy currents are induced in a conductor, when an electromagnetic 
field is incident upon it. Usually these phenomena result in flow of currents 
in a region near the boundary. A common way of treating these problems, 
in particular these induced currents, relies on estimating an appropriate 
skin depth on the conductor. Once the skin depth is determined, the as- 
sociated fields are obtained through an equivalence principle by imposing a 
boundary condition on the region of current formulation. Such a treatment 
may be found in [12]. The estimated skin depth (given at the end of the 
paper) is obtained by the penetration of fields into an infinite planar slab. 
However, if the conductor happens to have geometry other than that of a 
planar slab, the procedure mentioned above will not work well. Moreover, 
the incident waves need not be restricted to time-harmonic waves in general. 
Pulse sources and surge waves which come from thunder or circuit breakers 
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arise in practice. Therefore, one must have a procedure which reflects the 
transient behavior of eddy currents. 

Along this line, an effort was begun in [1], in which a frequency do- 
main problem was considered. This procedure presented some possibilities 
for further improvements and the feasibility of solving this class of prob- 
lems in three dimensions. The work in [3] dealt with continuation of [1] 
in conjunction with finite element methods. Unfortunately, only theoretical 
details are available in this work. Moreover, the treatment in the three- 
dimensional problem was restricted to time-harmonic fields. However, it 
was established in [4] and [5] that the finite difference methods are effective 
even for frequency domain problems, i. e. , pseudo-time marching. Our goal 
is to use finite difference methods to solve these problems, in particular, the 
two-dimensional problem in the time domain. In the engineering literature, 
for example, in [6], the author pointed out the computational difficulties 
in the problem. To motivate the ideas, the problem under consideration is 
presented in the next section. 

2 Formulation of the Problem 

Let fi be the cross-section region in the x-y plane of a conducting cylin- 
der with finite conductivity cr m , permeability n m , and permittivity e m . The 
generators of the cylinder are parallel to the z axis. The ambient medium 
is air with permeability fi a and permittivity e a . The conductivity of the air 
is neglected because it does not play a major role in this type of calcula- 
tions. If desired, it could be incorporated into this model, in which case the 
following treatment would still be valid. The incident field has the form 

E, = Ei(x,y,t) k, (1) 

H,- = Hn(x,y,t)i + H i2 (x,y,t)j . (2) 

A typical situiation is depicted in Fig. 1. The electric field is parallel to the 
axis of the cylinder and the magnetic field is transverse (TM) and in the 
x-y plane. Writing Maxwell’s equations, assuming all material properties 
are constant and isotropic, one obtains two sets of equations that are valid 
inside and outside the conductor respectively. 

In SIa (Air) 

dHi dH2 
dx ^ dy 


( 3 ) 
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dE_. _ he. 

dy 1 dx^ 
dH 1 _8H i 
dx dy 

In ft (Conductor) 

dHi dll 2 
dx dy 
dE_, _ dE_. 
dy 1 dx^ 
dHi dHi 
dx dy 


f d^. dH 2 . \ 

“*Ur 1 + T aT J J ’ 

(4) 

dE 

€a dt • 

(5) 

o , 

(6) 

(dHi. , dHi.\ 

M at 1+ dt J j ’ 

(7) 

„ dE 

o m E + e m — . 

(8) 


The first term on the right side of (8) is the conduction current and the sec- 
ond term is the displacement current. In addition, the following conditions 
are needed to make the problem well-posed. 

On T (Boundary of ft) Tangential components of E and H are contin- 
uous, i.e., no surface current density is assumed. 


At infinity the scattered E and H both decay to zero (Radiation Condi- 
tion). 

The zero divergence and two-dimensionality of H in both regions suggest 
the existence of a scalar function $ such that 


H i 


H 2 


d<$ 


dy 

dy 


dx 


(9) 

( 10 ) 


Substituting these in equations (4) and (7), we see that 


E = 



( 11 ) 


The function is identified as the magnetic vector potential. In turn, 

using this scalar function, the problem can be rewritten as follows. 



The governing equations are 


d 2 v 

d 2 v 

d 2 y 

dx 2 dy 2 

- €aHa dt2 

d 2 V 

d 2 V 


dx 2 dy 2 



d 2 y 
1 at 2 


The interface conditions on T are 


Ma 


m- 

(£)'- 


Mr 


( 


/0£\ + 

UJ ’ 


££\ + 


in fly i, (12) 
in fl. (13) 


(14) 

(15) 


where the minus sign indicates transition from the exterior to the boundary 
and the plus sign indicates transition from the boundary to the interior O. 
Since the problem is posed in the time domain, we require appropriate initial 
conditions. 


The initial conditions are 


V(x,y, 0) = $,(x,y,0) , 

d$(x,y,0) _ 6>$,(x,y,0) 

dt dt 


Sommerfeld’s radiation condition is 


lim r 1 / 2 



= 0 , 


(16) 

(17) 


(18) 


where ^ s represents the scattered wave. 


3 Difficulties in the Problem 

The problem presented in equations (12) through (18) yields an in- 
terface problem with T being the interface. This model is also applicable 
in calculating wave scattering by lossy dielectrics and the nondestructive 
evaluation(NDE) of various kinds of composite materials. In the event of 
composite materials, n m and e m will no longer be constants. Rather, they 
will be functions of space. The frequency of the incident wave may range 
from low frequencies to microwave frequencies. 

In solving this problem, one encounters certain difficulties. They are 
highlighted as follows: 



1. Implementation of the radiation condition. 

2. Implementation of the interface conditions. 

3. The large number of grid points that is needed to capture the skin 
depth in the high-frequency case. 

4. Appropriate local boundary condition for the low-frequency cases, 
where the wavelengths are much longer than the diameter of the con- 
ductor. 

To further address and answer the difficulties, we require scaling of the 
problem. Our scaling is performed as follows. 

1. Scale $ by ^o, i. e. , $ = where is the amplitude of the 

incident wave 

2. Scale x and y by L, i. e. , x = x/T, y = y/i, where L is a length of 
the order of the scatterer’s diameter. 

3. Scale t by 1/w, i. e. , ? = ut, where u is the angular velocity of the 
incident wave. 


As a result, the problem in these new coordinates has the following form: 
Governing Equations 


d 2 ^ d 2r2 d 2 ^ 

^=2+ -ZZT = “ L Vata—. 


Ox 2 dy 2 


dt 


2 > 


dH d 2 v r2 2r2 d 2 ^ 

dx 2 + dt ~ uL Mm<Tm dt +uL * lm€m dt 2 


(19) 

( 20 ) 


Interface Conditions 


Initial Conditions 


Ha* = 


(21) 

II 

d¥* 

dn 

(22) 

$(x,i/,0) = 

*K*,y,o) , 

(23) 

oW 

= 

dW], 

- W fx ' v ' 0) ■ 

(24) 
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Radiation Condition 


lim r 1 / 2 




= 0 , 


(25) 


where k a = coLy/ii a € a . 

The difficulty with these problems can be simply explained now. For 
instance, if we consider a 60 Hz incident field and a conductor with a 
diameter of .3141593 m , the parameters that appear in equations (19) 
and (20) have the following values: k a = uLy/ti a e a = 2.513274 x 10 -7 , 
k' m = uLy/iimCyn = 2.513274 X 10~ 7 , = Lyjufi m e m = 26.55045, where 

we have used: /i a = /i m = 1.256637 X 10“ 6 , € a = e m = 8.841941 X 10 -12 , 
a m - 3.72 x 10 7 , and L = 0.2 in the MKS unit system. Thus, reconsid- 
ering equation (19), we see that it behaves like Laplace’s equation, while 
equation (20) behaves like a parabolic equation. These observations were 
made in [13]. This only yields an approximation to the model given in equa- 
tions (19) through (25). We believe that even though the value of k a is 
extremely small, the problem is still governed by a wave phenomenon. This 
is further justified by the results that we present later. Thus, we do not 
approximate the model, but rather present procedures to solve the problem 
as we have derived. This yields several difficulties, particularly in terms of 
the computation. Our treatment is presented in the coming sections. 

Thus retaining the wave nature in the exterior, we further scale the 
problem for the purpose of treatment of boundary conditions. We scale the 
time by k a . The result yield the following nondimensional problem that we 
numerically solve: 


Governing Equations 


where 


d 2 % d 2 V 

d 2 v 

(26) 

dx 2 dy 2 

~ df 2 ' 

d 2 V d 2 V 

2 M 2 d 2 V 

(27) 

dx 2 dy 2 

~ l m Qt + k m g t 2 > 

l 2 

‘m 

£ 

b 

£ 

CN 

3 

1 


k ’ 

rta 


k 2 

UJ 2 L 2 fi m C m 



" *2 ' 



and 



Interface Conditions 


= /*m* + , 

dv~ _ dv+ 

dn dn 


Initial Conditions 


= $,(x,y,0) , 

dV , . d^i, 

W ix ' y ’ 0) = ~af( x ’ y ’°) ■ 


Radiation Condition 


lim r 1 / 2 

r— +OO 


( d*, 
V dr 


+ 



= 0 , 


Note that the bar notation has been dropped. 


(28) 

(29) 


(30) 

(31) 


(32) 


4 Radiation Condition 


The radiation condition as given in equation (32) is hard to implement 
numerically. It must be approximated in a suitable way before it is imposed 
on a computation domain. For the waves with moderate to high frequencies, 
the forms of this condition are well known. A summary may be found in [7]. 
The procedure here is to write the far-held scattered wave in the form: 




f(k a (r - t )) 

VT 




(33) 


This expression is the extended D’Alembert’s principle and indicates the 
fact that all the scattered fields are outgoing. From this, as in [7], it is easily 
verified that satisfies 


d<S s dV,, 1 T 

— - 4 - H = O 

dr + dt + 2r 



(34) 


However, for low-frequency waves, equation (34) does not provide satis- 
factory solutions. This is due to the fact that the wavelengths are very long 
and the substantial distances for the far-field boundaries require at least 
one wavelength. This is not practical for numerical implementation. Thus 



we introduce a new concept of formulating near-field conditions. Again the 
procedure is analogous to the one used in (33) and (34), but equation (33) 
will be replaced by an outgoing near- field solution. 

As suggested in [2], for low-frequency waves, the solution is dominated 
by logarithmic terms in both space and wave number. Thus one can seek 
solutions of the form: 



« (a + b In r)f(k a (r - t )) , 

(35) 

where a, b 

are constants. From this one obtains 




(36) 

where 

A < r >"r( a + ilnr) ' 

(37) 

For low-frequency time-harmonic waves, 



a « b ^7 + In , 

(38) 

where 7 is 

Euler’s constant. Substituting in A(r), we obtain 



-^( r ) r [ -y — In 2 -f \nk a + lnr. 

(39) 


It is easy to verify that A(r) is positive in the low-frequency case. 

With equations (34) and (36) imposed on the exterior boundary, the 
problem becomes well-posed. The proof for the uniqueness of the solution 
is shown in the appendix A. 

5 Finite Difference Formulation of the Problem 

The finite difference method has been found effective in time dependent 
hyperbolic problems. In addition, the radiation conditions in the previous 
section can easily be integrated in the scheme. Finite difference formulae 
used in solving this problem are summarized as follows. Note that the spatial 
differences in both x and y directions are chosen to be equal, i. e. , Ax = Ay. 
The superscript *i’ represents the time coordinate and the subscripts j and 
‘k’ represent the grid coordinates in the x and y directions respectively. 
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In the exterior region, the following central-difference formula is used to 
discretize equation (26): 

- « j? + (fj)* + *h+. + *i*-i - 4 * h) • 

(40) 

In the interior region, the first-order time derivative d^/dt is imple- 
mented by a forward-difference method. The second-order derivatives are 
implemented by central-difference methods. As the result, the following 
difference formula is obtained for Eq. (27): 

* j? = «*b - »j? + + «U* + *W« + *h-. - «!*)• («) 

where 


a = 


b = 


c — 


2 + jfrAt 

K m 

1 + &-At 


K m 

1 ( At \2 

j^e! 

1 + $-A t 


On the radiation boundary, the total wave consists of two components: 
the incident wave $,■ and the scattered wave i. e. , 


® = Vi + y s - 


(42) 


For the far field, the scattered wave propagates in a direction which is very 
close to the radial direction. The origin is located at the center of the 
scatterer as shown in Fig. 2. Using the derivatives in Cartesian coordinates, 
the radiation condition has the following form: 


89, 

dx 


COS0 + 


dy 


. a , 89, 9, 

s,n0 + — + 27 = 0 - 


(43) 


Finite difference methods usually yield difficulties at the corners of the 
rectangular radiation boundary. These are due to conflicts between the dif- 
ference formulae at the corners. As a result, these corners become sources of 
instability. To avoid corner problems, Engquist and Majda in [8] proposed 



a different boundary condition be used for the nearest two grid points to 
the corners. Unfortunately in this way, the implementation of the radia- 
tion boundary condition becomes more complicated. In our experiment, a 
method using a smooth transition was developed to solve corner problems. 
The approach is depicted in Fig. 3 and explained as follows. On the bound- 
ary A\ to A 2 , d^ s /dx is implemented by a forward-difference formula and 
d^a/dy is implemented by a backward-difference formula. On the boundary 
Bi to B 2 , dVa/dx and d^ s /dy are both implemented by backward-difference 
methods. It seems that implementation of dV 3 /dx on two boundaries will 
conflict at the intersection grid C. But since the term is multiplied by cos6 
and its value is 0 at C, the conflict can not take effect and corner problems 
are avoided. 

The difference formulae on the radiation boundary are summarized as 
follows. 


At Aq 


On A\ — A 2 


At C 


On B\ — B 2 


(*.)&' = «(*.)U-*t(«.)S + u -<*.&] < 44 ) 

*5? = (*!)#+ < 45 ) 

(*.);? = <■(*.)!* - - («.)hl 

-cl(*4» - ( 4 °) 

*;i' = mii' + (*.)&* w 

= (*•),¥ + (*•)&* < 40 > 

(«.)&' = «(*.)!» - *[(«.)$* - 

- («.)h-ll (“) 

= (*)&' + (».)£' < 51 ) 



(52) 

(53) 


At Bq 

(*.)&* = 

* 55 ,’ = <».•)&' + (* 4 ? 


Conditions on the lower boundary are handled by symmetry as the upper 
boundary. The coefficients a, b, and c are obtained in the following way: 

At 

2rAx ’ 

At 

b = cosd - — , 

Ax 

c = sinO^- . 

Ax 

The implementation of equation (36) is similar to the above. 

In the neighborhood of the interface, numerical implementation of the 
interface conditions is complicated. In reference [11], interface conditions are 
implemented by a cell integration method. Here we will propose a simple 
way to solve it directly using the finite difference method. 

First, for an interface boundary with some extent of curvature, we can 
approximate it by a polygon with its vertexes occupying the regular grids. 
For example, the circle T in Fig. 4 can be approximated by the polygon T . 
Then the field computation can be classified into two categories. 

1. For grids in the exterior region, the wave equation in the air is used. 

2. For grids in the interior region, the wave equation inside the conductor 
is used. 


For grids on the boundary, to obtain an accurate result, the following are 
suggested: 

1. When the grid separation is much larger than the estimated skin depth, 
the interior wave equation is used. 

2. When the grid separation is equal to or smaller than the estimated 
skin depth, the exterior wave equation is used. 

Secondly, the condition 

dv+ 

dn dn 


( 54 ) 



can be transformed into: 


BV- „ , 89- . a 89+ a , B9+ . 

— — cosd H — - — sin# = — — cost) -\ — - — stnd , 
ox dy ox dy 


(55) 


where 6 is the angle between the normal direction and the x axis. Since 
cosd and sinO are linearly independent, we have the following relations on 
the interface boundary: 


89~ _ 89+ 

dx dx ’ 

89~ _ 89+ 

dy dy 

With these, when the field computation involves a grid in the other region, an 
equivalent value at that grid can be obtained for computation. For example 
, the field computation at the grid A in Fig. 5 belongs to the category of 
the air. But a grid B in the conductor is involved in the difference formula 
for A as shown below: 




n +1 - 


= 2 n - ¥ j- 1 + (^) + 9k + - 4$ji) , 


(58) 


where ^ B is the equivalent value of $ at B. It can be obtained in the 
following way: 

+ (! - • ( 59 ) 

Substituting into the equation above, we obtain 

((H) i +9 i c+9 i D +9 i E -(3+l/fi r )9 i A ) . (60) 

Similarly, using the first interface condition, the field at the grid 1 can be 
calculated from the following: 

9\+ l = a9[ - 6* j- 1 + c((9iy/nr + (?$ y/Hr + n + n- 4*i) , (61) 

where a, b, and c are the same constants as those in equation (41). 

When the relative permeability fi r of the conductor is 1, approximating 
the interface boundary by some other geometry is not necessary. The field 
computation is classified in the same way as above. When a field computa- 
tion involves a grid in the other set, the field at that grid will be directly 



used for calculation because it equals to the equivalent value. For example, 
in Fig. 6, from the interface conditions, 

*5 = *B > ( 62 ) 

*5 - *i = (63) 

Thus *$£ = S&Q. As a result, the field at A can be obtained from the 
following formula: 

= 2n - ^ + (|^) 2 (V c + V D + ^ + * j. - 4 ’5'^) . (64) 

Similarly, the field at the grid 1 can be calculated from the following: 

V* 1 = a$j - by i" 1 + c(% + + % + % - 4*i) - (65) 

6 Numerical Experiments 

The numerical methods presented above were implemented in a com- 
puter program running on a STJN4/260 workstation. The conductor is a 
circular cylinder and is assumed to be aluminum. The grid coordinates with 
useful parameters are shown in Fig. 7. The incident wave is assumed to be 
a sinusoidal TM plane wave propagating in the x direction. The wave is as- 
sumed to be incident on the left boundary starting at t = 0. Therefore, the 
incident component at the grid point (j,k) in the application of the radiation 
condition is 

(^« )),k = cos[k a (iAt - jAx)] , (66) 

where k a = uLy/Ji^Z. The reason for this is that the incident wave which 
arrives at the point (j,k) is the incident wave at the point (0,k), j Ax time 
units earlier. 

The high frequency used in the numerical experiment is 2.387324 x 10 8 
Hz. With this frequency, parameters have the following values: 

• k a = 1 

• l 2 m = 2.804814 X 10 9 

• = 1 - 


13 



The numerical solution shown in Fig. 8 was obtained after ten periods from 
the beginning. The CPU time is 74 seconds. From the history plot at a 
point (67,34)(scc Fig. 13), the solution is found to be at steady state. 

The low frequency used in the experiment is 60 Hz. With this frequency, 
parameters have the following values: 

• k a = 2.513274 X 10 -7 

• 4 = 2.804814 X 10 9 

• k m = L 

The numerical solution using the near-field condition(see Fig. 10), is ob- 
tained after 2.513274 X 10 -4 period. The CPU time is about 1.7 hrs. From 
the history plot at a point (67,34)(scc Fig. 14), the solution in the exterior 
region is found to be in transient state but with a value almost approach- 
ing that at steady state. For reasons of comparison, the numerical solution 
using the far-field condition is also shown in Fig. 11. 


7 Comparison of Numerical and Analytic Solu- 
tion 


Since aluminum is a good conductor, In the high-frequency case, con- 
ductors such as aluminum behave close to a perfect conductor ( a m = oo). 
The analytic solution for the perfect-conductor problem is found in the lit- 
erature, [9] and [10], as follows: 




CO 

y: i n £ n cos(n(j>) 

71—0 


j n(^a^ ) 


0 , 


*?n (^a a ) 

n ( n\k a a) 



(67) 

( 68 ) 


where J n is the Bessel function of first kind and order n, II ^ is the Ilankel 
function of first kind and order n. £o = 1 an d e n = 2 forn > 1. This solution 
assumes a circular cylinder of radius a. With respect to these coordinates, 
the incident wave on the left boundary 


VjjK _ g— ik a XL g— ik a t f 


(69) 
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where xl is the distance from the center of the cylinder to the left boundary. 
In our coordinates, the incident wave on the left boundary is: 

$,• = e ~ ikat . (70) 


Therefore, the analytic solution must be modified as 


® (r,<M) 


V + (r,(f>, t ) 


OO 

T. i n e n cos(n4 > ) 

71=0 

e ik a XL e - ik at 


Jn(k a r ) 


0 . 


H { n ] (k a a) 



(71) 

(72) 


For t = 10 periods, the analytic solution (taking the first ten terms in the 
series) for the perfect-conductor case is plotted in Fig. 9. Comparing Fig. 9 
to Fig. 8, we find that the magnetic field distributions of the two solutions in 
air are in good agreement. The numerical solution also demonstrates a ten- 
dency of penetration of the wave into a good conductor. Deeper penetration 
occurs on the illuminated side than on the shadowed side. The difference 
between the two solutions near the artificial boundary comes from the fact 
that its distance from the center is not infinite. This trade-off, however, 
makes the numerical implementation on a computer possible and efficient. 

In the low-frequency case, the waves in both regions have the following 
analytic forms: 


OO 



= y i n £ n cos(ri(f))[J n (k a r ) + a n H^\k a r)] 

n=0 

gikaXL e ~ik at 

OO 

(73) 


= yi n e n cos(u(j))b n J n (pr)e lkaXL e~ ,kai , 

71=0 

(74) 


where p « Using interface conditions, the coefficients a n and 

b n are determined from the following equations: 

Jn(kaQ) *h = Mr^n*Ai(P®) t (75) 

K[j' n (Ka) + a n (H^)' (k a a)] = b n pj' n (pa) , (76) 


for n = 0,1,2, The first three coefficients in these two equations are 

calculated to be: 


a 0 


-1.027247 X 10 -2 - t‘9.971206 X 10 -2 



ai = -2.004432 X 10" 15 - i2.852592 X 10“ 14 
a 2 = -1.817832 X 10 -29 - il.288663 X 10 -28 
b 0 = -1.106546 X 10" 8 + i8.582522 X 10~ 9 
bi = 2.314644 X 10“ 14 + *3.575007 X 10 -14 
b 2 = 3.629113 X 10 -21 - *2.260305 X 10 -21 . 

Taking the first three terms in the series, we plot the analytic solution for the 
low frequency case in -Fig. 12. Comparing this to Fig. 11, we find that the 
numerical solution using the far-field condition has a large numeric error 
although its field pattern is similar to that of the analytic solution. The 
error comes from the position where the radiation condition is applied. This 
condition requires k a r to be large. But in our low-frequency experiment, k a r 
is very small (O(10 -6 )). As a result, the error appeares to be large. On the 
other hand, if we expand the computation domain to get a large k a r, the 
grid size will become too large to be implemented on a computer. Therefore 
a boundary condition which is appropriate to the case of low k a r is needed 
to solve the problem. 

Comparison of figures 10 and 12 shows that the numerical solution us- 
ing the near-field condition matches the analytic solution very well in the 
exterior region. The root mean square error for the whole computation do- 
main is calculated to be 1.772813 x 10 -3 . Both solutions at the points with 
r = 4,20,40,60, and 6 = 0, 7r/4, tt/2, 37t/4, tt are listed on Table 1. It shows 
that the fields on the illuminated and shadowed sides are nearly the same 
because of a strong diffraction for the low-frequency wave. The difference 
of the inside field between the two solutions is due to the fact that the nu- 
merical solution is in transient state while the analytic solution is at steady 
state. 

Vector potential is related to the scalar function by the following: 

A = k . (77) 

Using this relation, we obtain three-dimensional plots of vector potentials in 
various cases (see Fig. 15 through 18). In each plot, we see the magnitude of 
'3? versus x and y axes. From these plots, we can observe the variation of '5/ 
over the whole computational domain. Fig. 15 and Fig. 16 are the numerical 
and analytic solutions for the high-frequency case respectively. Fig. 17 and 
Fig. 18 are the numerical and analytic solutions for the low-frequency case 
respectively. We need these plots to check the numerical results near the 
interface boundary. According to the interface conditions, the fields are 



continuous on the boundary. Therefore must change smoothly near the 
boundary. No sharp jumps are allowed. Usually this is a strict test for a 
numerical result. Fig. 15 and Fig. 17 show that our numerical results pass 
the test. The values change smoothly and no jumps happen near the 
boundary. 

To prove that our implementation for an interface problem is correct, 
we need obtain a steady-state solution for the low-frequency case. Also to 
capture an accurate eddy current phenomenon, the grid separation must be 
chosen smaller than the estimated skin depth. This skin depth is defined in 
the case of a plane scatterer of infinite depth as: 


6 = 


Tt f UrnGn 


(78) 


where f is the frequency of the incident wave, and cr m are the permeability 
and conductivity of the scatterer respectively. In the computations that 
follow, we assume the conductor to be graphite with a m — 4x 10 4 5/Af(MKS 
units), /i m = /z a and e m = e a . The frequency of the incident wave is 6000 
Hz. In this case, the estimated skin depth is 3.248737 cm. The grid size 
is 25 x 25 with the diameter of the conductor occupying 21 grids. A grid 
separation of .6544985 cm is chosen. Therefore the grid separation is about 
1/5 of the estimated skin depth. The numerical results were obtained after 
one period from the beginning. It took 20 minutes CPU time on a CRAY- 
XMP-24. Fig. 19 and Fig. 20 show the results after one period of time. 
Comparing the corresponding analytic solutions in Fig. 21 and Fig. 22, 
the numerical results were found to be satisfactory. A clear example of eddy 
current phenomenon is demonstrated. Fig. 23 is the history at a point P(3,0) 
inside the conductor. A wave phenomenon is clearly observed in this figure. 
This justifies the comment we made previously. Both regions are governed 
by wave phenomena. Fig. 24 shows the numerical result for a quarter period. 
In this plot, a mesh of 61 grids by 61 grids, with the conductor’s diameter 
occupying 41 grids, is used. A grid separation of .3272493cm is chosen. A 
conductor with a m = 10 4 S/M, fi m = fi a and e m = t a is assumed. The 
frequency of the incident wave is also 6000 Hz. At a quarter period, the 
incident magnetic field becomes identically zero. However, from Fig. 24, we 
find a weak magnetic field still exists in the exterior region. At the beginning, 
the field penetrates from the exterior region into the interior region. During 
a quarter period, the interior field is so much stronger than the exterior 
field and thus it penetrates back to the exterior region through the interface 



conditions. Therefore the maximum value of W occurs inside the conductor 
at this time. For a general transient signal, the problem needs to be solved 
in the time domain so that this phenomenon can be understood further. 
Our procedure can deal with such situations. Corresponding calculations 
will be reported elsewhere . 

8 Concluding Remarks 

From our investigation, the following conclusions are obtained: 

1. Implementation of the interface problem in this paper is proved to 
be satisfactory for both high- and low-frequency cases, by numerical 
experiments. 

2. Implementation of the radiation condition is found effective in handling 
corner problems. 

3. To solve the low-frequency eddy current problem, a boundary condi- 
tion which is appropriate to the case of low k a r is needed. The bound- 
ary condition obtained has been shown to work very well in numerical 
experiments. 

4. To obtain an accurate eddy current representation, the grid separation 
near the interface should be smaller than the estimated skin depth in 
(78). 

5. In the low-frequency eddy current problem, our finite difference method 
calls for extremely fine time steps. As a result for good conductors 
such as aluminum, calculation of penetration of the field takes rela- 
tively larger time than that for conductors such as graphite. Thus an 
acceleration scheme is preferred so that calculations can be performed 
inexpensively. 



A Appendix: Uniqueness of the Solution 

The problem considered can be written in the following form: 


n 

= 

in i, 

(79) 


= A$ 

in f 

(80) 


= $ 

on r, 

(81) 

d\$s dxjfi 

d 9 


(82) 

dn ^ dn 

dn 

on r, 

■$ s (x,y,0),y s t (x,y,0) 

= 0 

in f Ia, 

(83) 

V(x,y,0),y t (x,y,0) 

= 0 

in f 2, 

(84) 

+ A(r)$ s 

= 0 

on r e . 

(85) 


Here A is the Laplacian operator. Subscripts l t\ ‘tt’, and V represent 
partial derivatives. Superscripts ‘s’ and T represent scattered wave and 
incident wave respectively. T e is the radiation boundary as shown in Fig. 
25. Here we assume F e to be a circle. A(r) = 1/2 r for the high-frequency 
case and 



1 * 

7 — ln2 + hxk a + lnr. 


( 86 ) 


for the low-frequency situation. 

First, assume that the problem has two solutions U\(x,y,t ) and U 2 (x,y,t ) 
and v(x , y, t ) = U 2 (x,y, t ) — ui(x, y, t ). Then 


Vtt 

= Av 

in 

(87) 

+ l 2 m v t 

= Av 

in $1, 

(88) 

v~ 

= v + 

on r, 

(89) 

dv~ 

dn 

dv + 

dn 

on T, 

(90) 

v(x,?/,0),u t (x,y,0) 

= 0 

in SIa, Q, 

(91) 

v r + V t + A(r)v 

= 0 

on r e . 

(92) 


We construct a function 

(v 2 + v« • \7v)dn A + 1 J J Q ( k ™ v ? + • vv)dn , (93) 


where V is the gradient operator and • is the operator for the inner product. 
Taking the derivative of E(t) with respect to t, we obtain 



E\t) = / J u (»*»« + V» • V^t)^ + J J^(k 2 m v t v tt + \?v -\7v t )dn, 

= J (v t Av + • S7vt)dSl A + J + V*> • Vvt - lm v t) d & 

= J (V • v t V v) dfi A + J /jV ■ t>« V » - 


/ v t V v ’2I^+ / VtSyv-nds— / / l^v^dQ, 

Js A Jr J Jn 

Jr V t v n ds - ^ v f v„d.s + ^ u t v n ds - J J Q l m v tdft 

Jr v t v n ds - J J q ImVtdSl 

I v t v r ds — I I lltfdQ. 

J r e i Jo 

^ v t (-Wt - A(r)v)ds - J J Q l m v tdtt ■ 


Thus we have 


(94) 


E\t) + J A(r)vv t ds = — J v^ds — jjj^dSl. (95) 

The right side of the equation above is less than or equal to zero. Then 

E'(t ) + §iJ ^A(r)v 2 ds < 0 . (96) 


From the initial condition, u(x,t/, 0) = 0, t7*(x,y, 0) = 0, and E( 0) = 0. It 
follows that 

E{t) + / ]-A(r)v 2 ds < 0 , (97) 

J r e 4 

so that 

E(t ) < — J I A(r)v 2 ds . (98) 

From the structure of A(r) (section 4) for both high- and low-frequency 
waves, A(r) > 0 on T e . As a result, E(t) < 0. But from the definition, 
E(t) > 0. Thus the possibility is E(t ) = 0. This leads to vt(x,y,t) = 0, i.e., 
v is constant for the time. Since v(x, y , 0) = 0, v(x, y, t) = 0 for all time and 
all space. This means ui(x,y,t) = U 2 (x,y,t). Therefore we conclude that 
the problem has a unique solution. 
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Table 1 . Solution Comparison for Low Frequency 


Radius 

Grid 

Numerical 

Analytic 

20 

( 80, 0) 
( 74, 14) 
( 60, 20) 
( 45, 14) 
( 40, 0) 

0 .8035427D-01 
0 . 7969991D-01 
0 .8028407D-01 
0 . 8197310D-01 
0 .8021339D-01 

0 . 7879715D-01 
0 . 7815606D-01 
0 . 7879725D-01 
0 . 8042111D-01 
0 . 7879735D-01 

40 

(100, 0) 
( 88, 28) 
( 60, 40) 
( 31, 28) 
( 20, 0) 

0 . 1246082D+00 
0 . 1239659D+00 
0 . 1245336D+00 
0 .1250927D+00 
0 .1244591D+00 

0 . 1227901D+00 
0 . 1221490D+00 
0 . 1227901D+00 
0 . 1232822D+00 
0 . 1227902D+00 

60 

(120, 0) 
(102, 42) 
( 60, 60) 
( 17, 42) 
( 0, 0) 

0 .1503752D+00 
0 .1498472D+00 
0 .1503764D+00 
0 .1505956D+00 
0 . 1503784D+00 

0 . 1485242D+00 
0 . 1478831D+00 
0 . 1485243D+00 
0 . 1486387D+00 
0 . 1485243D+00 

i 

4 

( 64, 0) 

( 62, 2) 
( 60, 4) 

( 57, 2) 

( 56, 0) 

0 .1710767D-13 
0 .1975745D-18 
0.1719809D-13 
0.3046863D-14 
0 .1731365D-13 

0 . 2201674D-04 
-0. 112747 4D-0 5 
0 . 2201422D-04 
0 . 10 69381D-04 
0 . 2201171D-04 


*In Double Precision expression. 





FIGURE 1. - TYPICAL SITUATION OF THE PROBLEM. 


FIGURE 2. - COORDINATES OF THE COMPUTATIONAL DOMAIN. 


F.D.: FORWARD DIFFERENCE 
B.D.: BACKWARD DIFFERENCE 
N. : NO EFFECT 
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FIGURE 5. - GRID POINTS NEAR THE INTERFACE BOUNDARY WHEN M r *1. 
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FIGURE 6. - GRID POINTS NEAR THE INTERFACE BOUNDARY WHEN H r =1. 



(120,60) 


( 120 , 0 ) 


(120,-60) 



MAGNETIC FIELD LINE, NUMERICAL 


FIGURE 7. - GRID SIZE FOR THE NUMERICAL EXPERIMENT 
IN SECTION 6. 


FIGURE 8. - CONTOUR PLOT OF THE NUMERICAL SOLUTION FOR THE 
HIGH-FREQUENCY CASE, f = 2.387324 X 10° Hz; GRID SIZE AND 
OTHER PARAMETERS ARE SHOWN IN FIG. 7; RESULT WAS OBTAINED 
AFTER 10 PERIODS FROM THE BEGINNING. 
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'LOT OF THE NUMERICAL SOLUTION. USING FIGURE 12. - COI 

iRY CONDITTION, FOR THE LOW-FREQUENCY LUTION OF FIG 

!ID SIZE AND OTHER PARAMETERS ARE SHOWN 
iS OBTAINED AFTER 2.513274 X 10 ^ PER- 

;ng. 























MAGNETIC FIELD LINE, NUMERICAL 

FIGURE 24.- CONTOUR PLOT OF THE NUMERICAL SOLUTION AT A 
QUARTER PERIOD, f = 6000 Hz; 0 M = 10 +f| S/M; U M = M a ; 
E m = £ a ; GRID SIZE IS 61 x 61 WITH THE DIAMETER OF THE 
CONDUCTOR OCCUPYING 41 GRIDS; GRID SEPARATION IS 
.3272493 CM. 
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